Capuchins Movements Animations in plantations

GIFs segmenting path with Hidden Markov Models

Index AKDE HMM Daily Travels

1 R Packages

Here I am using different R packages to manage spatial data (sf), make the estela behind the animated trajectories, and gganimate to unified ggplots as a GID file.

#devtools::install_github("16EAGLE/moveVis")
library(sf)
library(ggspatial)
library(gganimate)
library(tidyverse)
library(dplyr)
library(ggforce)
library(zoo)

2 Processed data from HMM

Using Capuchins movement data, specifically I am using processed data obtained after decoding behavioral states using Hidden Markov Models. Using viterbi algorithm I created the most likely sequence of states for each location. The idea is to create an animation to reflect how behavioral states (Resident and Transit) overlap with habitat types within plantation landscapes. It is expected that Resident is linked to native forest fragments, as this movement states is linked to resting, feeding and stationary foraging behaviors. Transit, a directonal and faster movement will be common in pine stands, as capuchins transpase this habitat to connect between high-resource areas (fragments). Low and high data.frames include months in which pine bark-stripping is less and more intense, respectively.

low <- readRDS("low_anim.rds")
high <- readRDS("high_anim.rds")

We have to select the columns that will be useful for our propose. there are a lot of other column that only will obsctale and delay the creation of the animation.

low <- low %>% select (ID, lat, lon, t, step, angle, tipouso_grouped, state2)
high <- high %>% select (ID, lat, lon, t, step, angle, tipouso_grouped, state2)

We will start with Columba group. Artificially to make the trajectories smoother, we will ingnore nighttime gaps (collares were program to avoide nigthtime data collection, capuchins do not move during night)… as GIF will be use timestamps as the criteria for changing frames.

colu <- high[high$ID == unique(high$ID)[2], ]


n <- nrow(colu)

colu$t2 <- seq(
  from = as.POSIXct("2022-03-18 09:30:00", tz = "UTC"),
  by = "30 mins",
  length.out = n
)

Animation are heavy to run, so I will select only a few rows (250), representing approximately nine days of movement. I will create the Simple Feature object to then align the CRS with the environment shape file obtained from Timber companies.

colu <- colu[200:550, ]

colu_sf <- st_as_sf(colu,
                    coords = c("lon", "lat"),
                    crs = 4326) 

2.1 Adding Plantation landscape ShapeFile

Shape fine from Timber companies. Plantations are mostly constitued of Pinus taeda trees but there are other species as well, we will group all as plantations. BNSM is Bosque Nativo sin Manejo (Native Forest without managment).

shp <- st_read("UNION.shp") %>%
  st_transform(crs = st_crs(colu_sf)) %>% 
  mutate(tipouso_grouped = ifelse(tipouso == "BNSM", "BNSM", "plantations"))

## Reading layer `UNION' from data source 
##   `C:\Users\Usuario\Desktop\Doctorado\movement_anim\anim_movement\UNION.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 16315 features and 25 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 729964.5 ymin: 7069500 xmax: 799942.6 ymax: 7150257
## Projected CRS: WGS 84 / UTM zone 21S

To have a better look of the relevant area (like zoom) we will create a box basing on the limits imposed by the locacion of Columba (during the nine selected days in this case). I will add a margin 0.1 grades (we are in WGS84) to have some air.

bbox_colu <- st_bbox(colu_sf)

# Expand a little bit... 

margin <- 0.01


bbox_vals <- c(
  xmin = as.numeric(bbox_colu["xmin"]) - margin,
  ymin = as.numeric(bbox_colu["ymin"]) - margin,
  xmax = as.numeric(bbox_colu["xmax"]) + margin,
  ymax = as.numeric(bbox_colu["ymax"]) + margin)

# Now we will use the expanded box

bbox_expanded <- structure(
  bbox_vals,
  class = "bbox",
  crs = st_crs(colu_sf))

# we have to make it Simple feature

bbox_poly <- st_as_sfc(bbox_expanded)

The plantation shapefile (shp) covers a much larger area than we need. To improve rendering speed and clarity in the animation, we crop the shapefile to include only the bounding box defined around Columba’s movement data. This reduces complexity and ensures our animation focuses on the relevant region. Before cropping, we apply two important steps to avoid issues with geometry operations:

1. Disable the use of the s2 geometry engine. The s2 engine is stricter and may fail when performing certain operations (like cropping or intersections) if geometries are not perfectly valid or cross the dateline.
2. Make sure all geometries in the shapefile are valid. This avoids errors during cropping, especially if the shapefile has complex or self-intersecting polygons.

sf::sf_use_s2(FALSE)

shp <- st_make_valid(shp)
shp_crop <- st_crop(shp, bbox_poly)

Here I prepare the Columba dataset by removing the spatial geometry and adding the X and Y coordinates as regular columns. This will make it easier to manipulate the data for visualization and animation.

2.2 Data.frame processing

colu_df <- colu_sf %>%
  st_drop_geometry() %>%
  cbind(st_coordinates(colu_sf))

2.2.1 Normalize the state variable for smoother animation

We convert the state variable into numeric format (just in case it’s stored as a factor or character). Then we rescale it to a 0–1 range usingscales::rescale, which will allow us to smoothly interpolate behavioral states in the animation (e.g., create gradients between states for smoother transitions).

colu_df$state2 <- as.numeric(as.character(colu_df$state2))

colu_df$state2_cont <- scales::rescale(colu_df$state2, to = c(0, 1))

2.2.2 Segment the movement into steps for animation

For the animation, we segment the trajectory into steps between each consecutive location. For each point, we find the next point (lead), and store the X/Y end coordinates and time. We also assign a unique frame number (row index), which will help animate the movement step by step.

colu_df_seg <- colu_df %>%
  arrange(t2) %>%
  mutate(
    xend = lead(X),
    yend = lead(Y),
    t2_end = lead(t2),
    state2_end = lead(state2_cont),
    frame = row_number()
  ) %>%
  filter(!is.na(xend))

2.2.3 Interpolate between GPS fixes for smoother paths

Here we interpolate the GPS trajectory to create a smoother animation path. We fill in intermediate steps every 0.25 index between each original GPS fix. This is useful for creating a smooth “flow” in the animation rather than jumping between the actual observed fixes. We interpolate time, location, and behavioral state. Finally I didnt we this, but could be explored in the future.

interp_df <- colu_df %>%
  arrange(t2) %>%
  mutate(index = row_number()) %>%
  complete(index = seq(min(index), max(index), by = 0.25)) %>%
  arrange(index) %>%
  mutate(
    t2 = zoo::na.approx(t2),
    X = zoo::na.approx(X),
    Y = zoo::na.approx(Y),
    state2_cont = zoo::na.approx(state2_cont)
  ) %>%
  drop_na()

To create a dynamic trail that follows the animal as it moves, we extract the N most recent movement segments for each animation frame. These will form the “tail” behind the animal.n_trail defines how many past segments will be shown at each frame. The older the segment, the more it will fade in the visualization. For each frame (i.e., step in the movement), we look backward up to n_trail steps and collect the corresponding segments. We also calculate the “trail age”, which tells us how many frames ago each segment occurred. This is used to control alpha (transparency) so older segments fade out gradually.

Optional: To make the movement trail look smoother and more natural, we use Bézier curves. Each movement segment (from a start point to an end point) is turned into a curve by adding a control point in the middle. This avoids sharp, angular turns. Bézier curves need three points: - Start (x1, y1) - Control point (x2, y2): midpoint between start and end - End (x3, y3)

#install.packages("ggforce")

n_trail <- 12

trail_df <- purrr::map_dfr(unique(colu_df_seg$frame), function(f) {
  colu_df_seg %>%
    filter(frame <= f & frame > f - n_trail) %>%
    mutate(trail_age = f - frame, current_frame = f)
}) %>%
  mutate(state_label = factor(round(state2),
                              levels = c(1, 2),
                              labels = c("Resident", "Transit")))

bezier_df <- trail_df %>%
  mutate(
    ctrl_x = (X + xend) / 2,
    ctrl_y = (Y + yend) / 2
  ) %>%
  select(current_frame, trail_age, state_label, 
         x1 = X, y1 = Y,
         x2 = ctrl_x, y2 = ctrl_y,
         x3 = xend, y3 = yend) %>%
  pivot_longer(cols = c(x1, x2, x3), names_to = "x_type", values_to = "x") %>%
  pivot_longer(cols = c(y1, y2, y3), names_to = "y_type", values_to = "y") %>%
  filter(substr(x_type, 2, 2) == substr(y_type, 2, 2)) %>%
  group_by(current_frame) %>%
  mutate(order = row_number()) %>%
  filter(n() == 3) %>%  
  ungroup()

3 Animations

3.1 Columba group animations

This is the first animations, simple and without home ranges and core areas.

ggplot() +
  geom_sf(data = shp_crop,
          aes(fill = tipouso_grouped),
          color = NA,
          alpha = 0.8) +

geom_bezier(data = bezier_df,
            aes(x = x, y = y, group = current_frame, color = state_label,
                alpha = 1 - trail_age / n_trail),
            size = 4,
            lineend = "round") +

  geom_segment(data = trail_df,
               aes(x = X, y = Y,
                   xend = xend, yend = yend,
                   group = current_frame,
                   color = state_label,
                   alpha = 1 - trail_age / n_trail),
               size = 3,
               lineend = "round") +

  geom_point(data = trail_df %>% filter(trail_age == 0),
             aes(x = xend, y = yend, color = state_label),
             size = 4) +

  scale_fill_manual(values = c("BNSM" = "olivedrab3", "plantations" = "saddlebrown")) +
  scale_color_manual(values = c("Resident" = "red3", "Transit" = "midnightblue")) +
  guides(alpha = "none") +

  coord_sf(xlim = c(st_bbox(bbox_poly)["xmin"], st_bbox(bbox_poly)["xmax"]),
           ylim = c(st_bbox(bbox_poly)["ymin"], st_bbox(bbox_poly)["ymax"])) +

  theme_minimal(base_size = 28) +
  labs(fill = "Uso del suelo",
       color = "Estado HMM") +

  transition_manual(current_frame) +
  ease_aes("cubic-in-out")

To save it.

animate(
  last_plot(),
  fps = 20,
  duration = 10,
  width = 900,
  height = 650,
  renderer = gifski_renderer("mono_curvo_con_estela2.gif")
)

3.1.1 Adding background colors

Here I am adding manually some saquares to fill spaces.. like water or human settlments. Is not the best way but I am avoiding to modify the shapefile.

make_rect <- function(xmin, xmax, ymin, ymax, crs) {
  st_as_sfc(st_bbox(c(
    xmin = as.numeric(xmin),
    xmax = as.numeric(xmax),
    ymin = as.numeric(ymin),
    ymax = as.numeric(ymax)
  ), crs = crs))
}

bbox <- st_bbox(shp_crop)

xmid <- mean(as.numeric(c(bbox["xmin"], bbox["xmax"])))
ymid <- mean(as.numeric(c(bbox["ymin"], bbox["ymax"])))

panel_lb <- make_rect(bbox["xmin"], xmid, bbox["ymin"], ymid, st_crs(shp_crop)) 
panel_lt <- make_rect(bbox["xmin"], xmid, ymid, bbox["ymax"], st_crs(shp_crop))  
panel_rb <- make_rect(xmid, bbox["xmax"], bbox["ymin"], ymid, st_crs(shp_crop))  
panel_rt <- make_rect(xmid, bbox["xmax"], ymid, bbox["ymax"], st_crs(shp_crop))  

3.1.2 Adding Home ranges and Core areas

In this section, I add the shapefiles representing the monthly home ranges and core areas estimated using autocorrelated kernel density estimators (AKDE). These shapefiles were generated using a custom function called sim_akde(). Detailed explanation of how the AKDEs were computed—including model selection, movement fitting, and range estimation, is in this github.io page

These home range estimates are used to contextualize the movement trajectories, showing how behavioral states relate to space use across different habitat types. Note that movement represent only nine days and home range and core areas were calculated for whole month (March).

hr_shp <- st_read("Columba_2022_mar_akde_95.shp") %>%
  st_transform(crs = st_crs(shp_crop)) 

## Reading layer `Columba_2022_mar_akde_95' from data source 
##   `C:\Users\Usuario\Desktop\Doctorado\movement_anim\anim_movement\Columba_2022_mar_akde_95.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 748830.3 ymin: 7141517 xmax: 751306.7 ymax: 7143789
## Projected CRS: WGS 84 / UTM zone 21S

ca_shp <- st_read("columna_2022_mar_akde_50.shp") %>%
  st_transform(crs = st_crs(shp_crop)) 

## Reading layer `columna_2022_mar_akde_50' from data source 
##   `C:\Users\Usuario\Desktop\Doctorado\movement_anim\anim_movement\columna_2022_mar_akde_50.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -54.51543 ymin: -25.81328 xmax: -54.50594 ymax: -25.80697
## Geodetic CRS:  WGS 84

bbox <- st_bbox(bbox_poly)

x_label <- mean(c(bbox["xmin"], bbox["xmax"]))
y_label <- as.numeric(bbox["ymax"]) + 200 

3.1.3 First GIF

ggplot() +
  
    geom_sf(data = panel_lb, fill = "lightblue", color = NA) +
    geom_sf(data = panel_lt, fill = "grey12", color = NA) +
    geom_sf(data = panel_rb, fill = "lightblue", color = NA) +
    geom_sf(data = panel_rt, fill = "grey12", color = NA)+

  geom_sf(data = shp_crop,
          aes(fill = tipouso_grouped),
          color = NA,
          alpha = 1) +
  
  geom_sf(data = hr_shp, fill = NA, color = "gold", linewidth = 1.6) +
  geom_sf(data = ca_shp, fill = NA, color = "khaki", linewidth = 1.2, linetype = "dashed") +


  geom_bezier(data = bezier_df,
              aes(x = x, y = y, group = current_frame, color = state_label,
                  alpha = 1 - trail_age / n_trail),
              size = 5,
              lineend = "round") +

  geom_segment(data = trail_df,
               aes(x = X, y = Y, xend = xend, yend = yend,
                   group = current_frame,
                   color = state_label,
                   alpha = 1 - trail_age / n_trail),
               size = 5,
               lineend = "round") +

  geom_point(data = trail_df %>% filter(trail_age == 0),
             aes(x = xend, y = yend, color = state_label),
             size = 5.5) +
  
  annotate("label",
         x = x_label, y = y_label,
         label = "22 March 2022 to 31 March 2022\nGroup: Columba",
         hjust = 1.05, vjust = 1.5,
         size = 8,
         family = "sans",
         fontface = "bold",
         fill = "white",
         label.size = NA,  
         color = "black") +


  scale_fill_manual(
    values = c("BNSM" = "olivedrab3", "plantations" = "#7f5529"),
    labels = c("BNSM" = "Native Forest", "plantations" = "Pine stands"),
    name = "Habitat type",
    na.translate = FALSE) +

  scale_color_manual(
    values = c("Resident" = "midnightblue", "Transit" = "red3"),
    name = "HMM States") +

    guides(
    alpha = "none",
    linetype = guide_legend(
      override.aes = list(
      color = c("gold", "khaki"),
      size = 1.6))) +

  coord_sf(xlim = c(st_bbox(bbox_poly)["xmin"], st_bbox(bbox_poly)["xmax"]),
           ylim = c(st_bbox(bbox_poly)["ymin"], st_bbox(bbox_poly)["ymax"])) +

  theme_minimal(base_size = 28, base_family = "sans") +
  theme(
    axis.text.x = element_text(angle = 45, size = 21, color = "black", family = "sans"),
    axis.text.y = element_text(color = "black", size = 21, family = "sans"),
    axis.title.x = element_text(color = "black", size = 23, family = "sans"),
    axis.title.y = element_text(color = "black", size = 23, family = "sans"),
    text = element_text(color = "black", family = "sans"),
    legend.text = element_text(color = "black", size = 26, family = "sans"),
    legend.title = element_text(color = "black", size = 28, family = "sans")) +

  labs(fill = "Habitat type",
       color = "HMM States",
       x = "Longitude",
       y = "Latitude") + 
  annotation_scale(
  location = "bl",  
  width_hint = 0.4,
  bar_width = 1000,     
  unit = "m",           
  text_cex = 1.4,  
  pad_x = unit(0.6, "in"),
  pad_y = unit(0.6, "in")) +
  
  annotation_north_arrow(
    location = "tl", which_north = "true",
    style = north_arrow_fancy_orienteering,
    pad_x = unit(0.6, "in"),
    pad_y = unit(0.6, "in")) +

  transition_manual(current_frame) +
  ease_aes(y = 'bounce-in')

animate(
  last_plot(),
  fps = 30,
  duration = 16,
  width = 2800,
  height = 1450,
  res = 150,
  renderer = gifski_renderer("mono_curvo_con_estela_hr.gif"))

3.1.4 Second GIF

Adding caption and improving ranging areas visualization

ggplot() +

  geom_sf(data = panel_lb, fill = "lightblue", color = NA) +
  geom_sf(data = panel_lt, fill = "grey12", color = NA) +
  geom_sf(data = panel_rb, fill = "lightblue", color = NA) +
  geom_sf(data = panel_rt, fill = "grey12", color = NA) +

  geom_sf(data = shp_crop,
          aes(fill = tipouso_grouped),
          color = NA,
          alpha = 1) +

  geom_sf(data = ca_shp, aes(linetype = factor("Core area", levels = c("Core area", "Home range"))), 
        fill = NA, color = "goldenrod1", linewidth = 2.9) +

  geom_sf(data = hr_shp, aes(linetype = factor("Home range", levels = c("Core area", "Home range"))), 
        fill = NA, color = "yellow", linewidth = 2.9) +

  geom_bezier(data = bezier_df,
              aes(x = x, y = y, group = current_frame, color = state_label,
                  alpha = 1 - trail_age / n_trail),
              size = 8,
              lineend = "round") +

  geom_segment(data = trail_df,
               aes(x = X, y = Y, xend = xend, yend = yend,
                   group = current_frame,
                   color = state_label,
                   alpha = 1 - trail_age / n_trail),
               size = 8,
               lineend = "round") +

  geom_point(data = trail_df %>% filter(trail_age == 0),
             aes(x = xend, y = yend, color = state_label),
             size = 8.5) +

  scale_fill_manual(
    values = c("BNSM" = "olivedrab3", "plantations" = "#7f5529"),
    labels = c("BNSM" = "Native Forest", "plantations" = "Pine stands"),
    name = "Habitat type",
    na.translate = FALSE) +

  scale_color_manual(
    values = c("Resident" = "midnightblue", "Transit" = "red3"),
    name = "HMM States") +

  scale_linetype_manual(
    values = c("Home range" = "solid", "Core area" = "solid"),
    name = "Ranging areas") +

 guides(
  alpha = "none",
  linetype = guide_legend(
    override.aes = list(
      linetype = c("solid", "solid"),
      color = c("goldenrod1", "yellow"),
      size = c(1.9, 1.9)
    ),
    title = "Ranging areas" )) +

  coord_sf(xlim = c(st_bbox(bbox_poly)["xmin"], st_bbox(bbox_poly)["xmax"]),
           ylim = c(st_bbox(bbox_poly)["ymin"], st_bbox(bbox_poly)["ymax"])) +

  
  ggtitle("Columba group, 22 March to 31 March 2022") +

 
  theme_minimal(base_size = 28, base_family = "sans") +
  theme(
    plot.title = element_text(hjust = 0.5, size = 29.5, face = "bold"),
    axis.text.x = element_text(angle = 45, size = 22, color = "black"),
    axis.text.y = element_text(size = 22,  color = "black"),
    axis.title.x = element_text(size = 26, color = "black"),
    axis.title.y = element_text(size = 26, color = "black"),
    text = element_text(color = "black"),
    legend.text = element_text(size = 29, color = "black"),
    legend.title = element_text(size = 31, color = "black"),
    plot.margin = margin(t = 10, r = 10, b = 150, l = 10),
    plot.caption = element_text(
  size = 20, hjust = 0.2,vjust = -4.5, 
  color = "black", 
  family = "sans",
  margin = margin(t = 80))) +

  labs(
    fill = "Habitat type",
    color = "HMM States",
    x = "Longitude",
    y = "Latitude",
    caption = "Movement states of black capuchin monkeys (*Sapajus nigritus*) in pine plantations of northeastern Argentina.\nThe trajectory is segmented based on behavioral states inferred from Hidden Markov Models.\nMonthly home ranges and core areas were estimated using 95% and 50% autocorrelated kernel density\n estimators (AKDE), respectively. The animation illustrates how movement segmentation is associated with habitat type:\n*Resident* behavior is more common in native forest, while *Transit* behavior predominates in pine stands. \n Lightblue represent water") +

  annotation_scale(
    location = "bl",
    width_hint = 0.4,
    bar_width = 1000,
    unit = "m",
    text_cex = 1.4,
    pad_x = unit(0.6, "in"),
    pad_y = unit(0.6, "in") ) +

  annotation_north_arrow(
    location = "tl", which_north = "true",
    style = north_arrow_fancy_orienteering,
    pad_x = unit(0.6, "in"),
    pad_y = unit(0.6, "in")
  ) +


  transition_manual(current_frame) +
  ease_aes(y = "bounce-in")

animate(
  last_plot(),
  fps = 30,
  duration = 16,
  width = 2100,
  height = 1550,
  res = 100,
  renderer = gifski_renderer("mono_curvo_con_estela_hr_ca.gif")
)

3.2 Gualicho group animations

Now we will do the same but with another group. The groups is in the south so it include some pine plantations from PINDO company, this could be a problem, the company shapefile is differently structure compared to that of Arauco.

pan <- high[high$ID == unique(high$ID)[1], ]

n <- nrow(pan)

pan$t2 <- seq(
  from = as.POSIXct("2022-08-18 09:30:00", tz = "UTC"),
  by = "30 mins",
  length.out = n
)

pan <- pan[150:450, ]

pan_sf <- st_as_sf(pan,
                    coords = c("lon", "lat"),
                    crs = 4326) 

shp_pan <- st_read("UNION_GUALICHO.shp") %>%
  st_transform(crs = st_crs(pan_sf)) %>%  # Alinear CRS
  mutate(tipouso_grouped = ifelse(tipouso == "BNSM", "BNSM", "plantations"))

## Reading layer `UNION_GUALICHO' from data source 
##   `C:\Users\Usuario\Desktop\Doctorado\movement_anim\anim_movement\UNION_GUALICHO.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 16318 features and 25 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 727994.5 ymin: 7069500 xmax: 799942.6 ymax: 7150257
## Projected CRS: WGS 84 / UTM zone 21S

bbox_pan <- st_bbox(pan_sf)

margin <- 0.01

bbox_vals <- c(
  xmin = as.numeric(bbox_pan["xmin"]) - margin,
  ymin = as.numeric(bbox_pan["ymin"]) - margin,
  xmax = as.numeric(bbox_pan["xmax"]) + margin,
  ymax = as.numeric(bbox_pan["ymax"]) + margin
)

bbox_expanded <- structure(
  bbox_vals,
  class = "bbox",
  crs = st_crs(pan_sf)
)

bbox_poly <- st_as_sfc(bbox_expanded)

cropping to fit Tilo’s area

sf::sf_use_s2(FALSE)

shp_pan <- st_make_valid(shp_pan)
shp_crop_pan <- st_crop(shp_pan, bbox_poly)

pan_df <- pan_sf %>%
  st_drop_geometry() %>%
  cbind(st_coordinates(pan_sf))

pan_df$state2 <- as.numeric(as.character(pan_df$state2))

pan_df$state2_cont <- scales::rescale(pan_df$state2, to = c(0, 1))

pan_df_seg <- pan_df %>%
  arrange(t2) %>%
  mutate(
    xend = lead(X),
    yend = lead(Y),
    t2_end = lead(t2),
    state2_end = lead(state2_cont),
    frame = row_number()
  ) %>%
  filter(!is.na(xend))

interp_df_pan <- pan_df %>%
  arrange(t2) %>%
  mutate(index = row_number()) %>%
  complete(index = seq(min(index), max(index), by = 0.25)) %>%
  arrange(index) %>%
  mutate(
    t2 = zoo::na.approx(t2),
    X = zoo::na.approx(X),
    Y = zoo::na.approx(Y),
    state2_cont = zoo::na.approx(state2_cont)
  ) %>%
  drop_na()

n_trail <- 10

trail_df_pan <- purrr::map_dfr(unique(pan_df_seg$frame), function(f) {
  pan_df_seg %>%
    filter(frame <= f & frame > f - n_trail) %>%
    mutate(trail_age = f - frame, current_frame = f)
}) %>%
  mutate(state_label = factor(round(state2),
                              levels = c(1, 2),
                              labels = c("Resident", "Transit")))

bezier_df_pan <- trail_df_pan %>%
  mutate(
    ctrl_x = (X + xend) / 2,
    ctrl_y = (Y + yend) / 2
  ) %>%
  select(current_frame, trail_age, state_label, 
         x1 = X, y1 = Y,
         x2 = ctrl_x, y2 = ctrl_y,
         x3 = xend, y3 = yend) %>%
  pivot_longer(cols = c(x1, x2, x3), names_to = "x_type", values_to = "x") %>%
  pivot_longer(cols = c(y1, y2, y3), names_to = "y_type", values_to = "y") %>%
  filter(substr(x_type, 2, 2) == substr(y_type, 2, 2)) %>%
  group_by(current_frame) %>%
  mutate(order = row_number()) %>%
  filter(n() == 3) %>%  
  ungroup()

make_rect <- function(xmin, xmax, ymin, ymax, crs) {
  st_as_sfc(st_bbox(c(
    xmin = as.numeric(xmin),
    xmax = as.numeric(xmax),
    ymin = as.numeric(ymin),
    ymax = as.numeric(ymax)
  ), crs = crs))
}


bbox_pan <- st_bbox(shp_crop_pan)

xmid_pan <- mean(as.numeric(c(bbox_pan["xmin"], bbox_pan["xmax"])))
ymid_pan <- mean(as.numeric(c(bbox_pan["ymin"], bbox_pan["ymax"])))

panel_lb_pan <-
  make_rect(bbox_pan["xmin"], 
            xmid_pan, 
            bbox_pan["ymin"],
            ymid_pan, 
            st_crs(shp_crop_pan)) 

panel_lt_pan <- 
  make_rect(bbox_pan["xmin"],
            xmid_pan,
            ymid_pan, 
            bbox_pan["ymax"], 
            st_crs(shp_crop_pan))  

panel_rb_pan <-
  make_rect(xmid_pan, 
            bbox_pan["xmax"],
            bbox_pan["ymin"],
            ymid_pan, 
            st_crs(shp_crop_pan))  

panel_rt_pan <- 
  make_rect(xmid_pan,
            bbox_pan["xmax"],
            ymid_pan,
            bbox_pan["ymax"], 
            st_crs(shp_crop_pan))  

3.2.1 Home ranges and core areas

hr_shp_pan <- st_read("gualicho_2022_ago_akde_95.shp") %>%
  st_transform(crs = st_crs(shp_crop)) 

## Reading layer `gualicho_2022_ago_akde_95' from data source 
##   `C:\Users\Usuario\Desktop\Doctorado\movement_anim\anim_movement\gualicho_2022_ago_akde_95.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 732491.5 ymin: 7118941 xmax: 736624.5 ymax: 7122846
## Projected CRS: WGS 84 / UTM zone 21S

ca_shp_pan <- st_read("gualicho_2022_ago_akde_50.shp") %>%
  st_transform(crs = st_crs(shp_crop)) 

## Reading layer `gualicho_2022_ago_akde_50' from data source 
##   `C:\Users\Usuario\Desktop\Doctorado\movement_anim\anim_movement\gualicho_2022_ago_akde_50.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 732596.2 ymin: 7119936 xmax: 733793.7 ymax: 7122528
## Projected CRS: WGS 84 / UTM zone 21S

bezier_df_pan <- bezier_df_pan %>% filter(!is.na(current_frame))
trail_df_pan <- trail_df_pan %>% filter(!is.na(current_frame))

3.2.2 Animation

ggplot() +

  geom_sf(data = panel_lb_pan, fill = "lightblue", color = NA) +
  geom_sf(data = panel_lt_pan, fill = "lightblue", color = NA) +
  geom_sf(data = panel_rb_pan, fill = "grey74", color = NA) +
  geom_sf(data = panel_rt_pan, fill = "grey74", color = NA) +

  geom_sf(data = shp_crop_pan,
          aes(fill = tipouso_grouped),
          color = NA,
          alpha = 1) +

  geom_sf(data = ca_shp_pan, aes(linetype = factor("Core area", levels = c("Core area", "Home range"))), 
        fill = NA, color = "darkorange", linewidth = 2.9) +

  geom_sf(data = hr_shp_pan, aes(linetype = factor("Home range", levels = c("Core area", "Home range"))), 
        fill = NA, color = "yellow", linewidth = 2.9) +

  geom_bezier(data = bezier_df_pan,
              aes(x = x, y = y, group = current_frame, color = state_label,
                  alpha = 1 - trail_age / n_trail),
              size = 8,
              lineend = "round") +

  geom_segment(data = trail_df_pan,
               aes(x = X, y = Y, xend = xend, yend = yend,
                   group = current_frame,
                   color = state_label,
                   alpha = 1 - trail_age / n_trail),
               size = 8,
               lineend = "round") +

  geom_point(data = trail_df_pan %>% filter(trail_age == 0),
             aes(x = xend, y = yend, color = state_label),
             size = 8.5) +

  scale_fill_manual(
    values = c("BNSM" = "olivedrab3", "plantations" = "#7f5529", "human settlments" = "grey74"),
    labels = c("BNSM" = "Native Forest", "plantations" = "Pine stands"),
    name = "Habitat type",
    na.translate = FALSE) +

  scale_color_manual(
    values = c("Resident" = "midnightblue", "Transit" = "red3"),
    name = "HMM States") +

  scale_linetype_manual(
    values = c("Home range" = "solid", "Core area" = "solid"),
    name = "Ranging areas") +

 guides(
  alpha = "none",
  linetype = guide_legend(
    override.aes = list(
      linetype = c("solid", "solid"),
      color = c("darkorange", "yellow"),
      size = c(1.9, 1.9)),
    title = "Ranging areas")) +

  coord_sf(xlim = c(st_bbox(bbox_poly)["xmin"], st_bbox(bbox_poly)["xmax"]),
           ylim = c(st_bbox(bbox_poly)["ymin"], st_bbox(bbox_poly)["ymax"])) +

  
  ggtitle("Gualicho group, 18 August to 27 August 2022") +

  theme_minimal(base_size = 28, base_family = "sans") +
  theme(
    plot.title = element_text(hjust = 0.5, size = 29.5, face = "bold"),
    axis.text.x = element_text(angle = 45, size = 22, color = "black"),
    axis.text.y = element_text(size = 22,  color = "black"),
    axis.title.x = element_text(size = 26, color = "black"),
    axis.title.y = element_text(size = 26, color = "black"),
    text = element_text(color = "black"),
    legend.text = element_text(size = 29, color = "black"),
    legend.title = element_text(size = 31, color = "black"),
    plot.margin = margin(t = 10, r = 10, b = 150, l = 10),
    plot.caption = element_text(
  size = 20, hjust = 0.2,vjust = -4.5, 
  color = "black", 
  family = "sans",
  margin = margin(t = 80))) +

  labs(
    fill = "Habitat type",
    color = "HMM States",
    x = "Longitude",
    y = "Latitude",
    caption = "Movement states of black capuchin monkeys (*Sapajus nigritus*) in pine plantations of northeastern Argentina.\nThe trajectory is segmented based on behavioral states inferred from Hidden Markov Models.\nMonthly home ranges and core areas were estimated using 95% and 50% autocorrelated kernel density\n estimators (AKDE), respectively. The animation illustrates how movement segmentation is associated with habitat type:\n*Resident* behavior is more common in native forest, while *Transit* behavior predominates in pine stands.") +

  annotation_scale(
    location = "bl",
    width_hint = 0.4,
    bar_width = 1000,
    unit = "m",
    text_cex = 1.4,
    pad_x = unit(0.6, "in"),
    pad_y = unit(0.6, "in") ) +

  annotation_north_arrow(
    location = "tl", which_north = "true",
    style = north_arrow_fancy_orienteering,
    pad_x = unit(0.6, "in"),
    pad_y = unit(0.6, "in")
  ) +

  transition_manual(current_frame) +
  ease_aes(y = "bounce-in")

animate(
  last_plot(),
  fps = 30,
  duration = 16,
  width = 2100,
  height = 1550,
  res = 100,
  renderer = gifski_renderer("mono_curvo_con_estela_hr_ca_gua.gif"))

3.2.3 High quality renderization

Biological Conservation for instance requires 300 DPI in animations, and more than 15 fps.

ggplot() +

  geom_sf(data = panel_lb_pan, fill = "lightblue", color = NA) +
  geom_sf(data = panel_lt_pan, fill = "lightblue", color = NA) +
  geom_sf(data = panel_rb_pan, fill = "grey34", color = NA) +
  geom_sf(data = panel_rt_pan, fill = "grey34", color = NA) +

  geom_sf(data = shp_crop_pan,
          aes(fill = tipouso_grouped),
          color = NA,
          alpha = 1) +

  geom_sf(data = ca_shp_pan, aes(linetype = factor("Core area", levels = c("Core area", "Home range"))), 
        fill = NA, color = "darkorange", linewidth = 1.8) +
  geom_sf(data = hr_shp_pan, aes(linetype = factor("Home range", levels = c("Core area", "Home range"))), 
        fill = NA, color = "yellow", linewidth = 1.8) +

  geom_bezier(data = bezier_df_pan,
              aes(x = x, y = y, group = current_frame, color = state_label,
                  alpha = 1 - trail_age / n_trail),
              size = 5.5,
              lineend = "round") +

  geom_segment(data = trail_df_pan,
               aes(x = X, y = Y, xend = xend, yend = yend,
                   group = current_frame,
                   color = state_label,
                   alpha = 1 - trail_age / n_trail),
               size = 5.5,
               lineend = "round") +

  geom_point(data = trail_df_pan %>% filter(trail_age == 0),
             aes(x = xend, y = yend, color = state_label),
             size = 6) +

  scale_fill_manual(
    values = c("BNSM" = "olivedrab3", "plantations" = "#7f5529"),
    labels = c("BNSM" = "Native Forest", "plantations" = "Pine stands"),
    name = "Habitat type",
    na.translate = FALSE) +

  scale_color_manual(
    values = c("Resident" = "midnightblue", "Transit" = "red3"),
    name = "HMM States") +

  scale_linetype_manual(
    values = c("Home range" = "solid", "Core area" = "solid"),
    name = "Ranging areas") +

 guides(
  alpha = "none",
  linetype = guide_legend(
    override.aes = list(
      linetype = c("solid", "solid"),
      color = c("darkorange", "yellow"),
      size = c(1.8, 1.8)),
    title = "Ranging areas")) +

  coord_sf(xlim = c(st_bbox(bbox_poly)["xmin"], st_bbox(bbox_poly)["xmax"]),
           ylim = c(st_bbox(bbox_poly)["ymin"], st_bbox(bbox_poly)["ymax"])) +

  ggtitle("Gualicho group, 18 August to 27 August 2022") +

  theme_minimal(base_size = 20, base_family = "sans") +
  theme(
    plot.title = element_text(hjust = 0.5, size = 22, face = "bold"),
    axis.text.x = element_text(angle = 45, size = 15, color = "black"),
    axis.text.y = element_text(size = 15, color = "black"),
    axis.title.x = element_text(size = 18, color = "black"),
    axis.title.y = element_text(size = 18, color = "black"),
    text = element_text(color = "black"),
    legend.text = element_text(size = 20, color = "black"),
    legend.title = element_text(size = 22, color = "black"),
    plot.margin = margin(t = 10, r = 10, b = 180, l = 10),
    plot.caption = element_text(
      size = 15,
      hjust = 0.5,
      vjust = -4.5,
      color = "black",
      family = "sans",
      margin = margin(t = 80))) +

  labs(
    fill = "Habitat type",
    color = "HMM States",
    x = "Longitude",
    y = "Latitude",
    caption = "Movement states of black capuchin monkeys (Sapajus nigritus) in pine plantations of northeastern Argentina.\nThe trajectory is segmented based on behavioral states inferred from Hidden Markov Models.\nMonthly home ranges and core areas were estimated using 95% and 50% autocorrelated kernel density\nestimators (AKDE), respectively. The animation illustrates how movement segmentation is associated with habitat type:\nResident behavior is more common in native forest, while Transit behavior predominates in pine stands.\n Lightblue represent water and Grey human settlements.") +

  annotation_scale(
    location = "bl",
    width_hint = 0.4,
    bar_width = 1000,
    unit = "m",
    text_cex = 1.2,
    pad_x = unit(0.6, "in"),
    pad_y = unit(0.6, "in") ) +

  annotation_north_arrow(
    location = "tl", which_north = "true",
    style = north_arrow_fancy_orienteering,
    pad_x = unit(0.6, "in"),
    pad_y = unit(0.6, "in")) +

  transition_manual(current_frame) +
  ease_aes(y = "cubic-in-out")

animate(
  last_plot(),
  fps = 30,
  duration = 16,
  width = 6300,      # 3 × original width
  height = 4650,     # 3 × original height
  res = 300,
  renderer = gifski_renderer("gualicho_300ppi.gif")
)